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ABSTRACT 



Aims. We aim to investigate the abundances of light deuterium-bearing species such as HD, H2D + and D2H + in a gas-grain chemical 
model including an extensive description of deuterium and spin state chemistry, in physical conditions appropriate to the very centers 
of starless cores. 

Methods. We combine a gas-grain chemical model with radiative transfer calculations to simulate density and temperature structure 
in starless cores. The chemical model includes new reaction sets for both gas phase and grain surface chemistry, including deuterated 
forms of species with up to 4 atoms and the spin states of the light species H 2 , HJ and HJ and their deuterated forms. 
Results. We find that in the dense and cold environments attributed to the centers of starless cores, HD eventually depletes from the 
gas phase because deuterium is efficiently incorporated to grain-surface HDO, resulting in inefficient HD production on grains for 
advanced core ages. HD depletion has consequences not only on the abundances of e.g. H 2 D + and D2H + , whose production depends 
on the abundance of HD, but also on the spin state abundance ratios of the various light species, when compared with the complete 
depletion model where heavy elements do not influence the chemistry. 

Conclusions. While the eventual HD depletion leads to the disappearance of light deuterium-bearing species from the gas phase in 
a relatively short timescale at high density, we find that at late stages of core evolution the abundances of H 2 D + and D 2 H + increase 
toward the core edge and the disributions become extended. The HD depletion timescale increases if less oxygen is initially present 
in the gas phase, owing to chemical interaction between the gas and the dust predecing the starless core phase. Our results are greatly 
affected if H 2 is allowed to tunnel on grain surfaces, and therefore more experimental data not only on tunneling but also on the O + H 2 
surface reaction in particular is needed. 

Key words. ISM: abundances - ISM: clouds - ISM: molecules - astrochemistry - radiative transfer 



1. Introduction 

Starless cores, sites for potential low-mass star formation, are 
condensations of dense and cold gas. Owing to the high den- 
sity and low temperature attributed to these objects, it is ex- 
pected that species containing heavy elements eventually de- 
plete onto grain surfaces because of their relatively high bind- 
ing energie s. Indeed, depletion of several ch emical species such 
as CO (e.g.lWillacv et al.lll998r.lCaselli et alJll999h and CS (e.g. 
' iTafalla et alj|2002t iPon et al.ll2009l) has been observed toward 
these objects. 

For a long time (since the early 1970s), gas-phase chemi- 
cal models were the main means of studies of chemical evolu- 
tion in starless cores. However, even though it has been known 
for a long time that H2 must be formed on grain surfaces 



dGould & Salpeteilll963l:lHollenbach & Salpeteill 19711) and con- 
sequently that surface chemistry should play an important role in 
the chemical evolution of these objects, it is only fairly recently 
(during the last 10 years or so) that models including the interac- 
tion between gas phase and grain surface chemistry have really 
taken over as the main means of numerical studies of the chem- 
istry in starless cores (lAikawa et al.l 120031 120051: iGarrod et alJ 
120071) . 

One of the main goals of simulating chemistry in starless 
cores is not only to explain the observed properties (such as line 
emission radiation from various species) of these cores, but also 
to search for new tracer species for the varying physical con- 



ditions. Particularly in the very centers of starless cores where 
density is high, species containing heavy elements are expected 
to be almost totally depleted onto grain surfaces, and it has 
been suggested that the light deuterium-bearing species H2D + 
and D2FP could serve as the main trace rs of these conditions 
dWalmslev et al.ll2004l:lFlower et al.ll2004l) . after observations of 
the ground state rotational transition of ortho-H2D + revealed 
an unexpected stro ng line toward the prestellar core L1544 
(ICaselli et al.ll2003J). However, in the so-c alled complet e deple - 
tion model of lWalmslev et all d2004 and iFlower et all (120041) . 
in which no heavy elements are present in the gas phase, a de- 
scription of grain surface chemistry is not included and it has re- 
mained somewhat unclear whether the abundances of H2D + and 
D2H + would be affected if one were to consider surface chem- 
istry as well. 

In this paper, we address this issue by studying the chem- 
istry in physical conditions appropriate to the centers of starless 
cores using a gas-grain chemical model. To this end, we have 
constructed new chemical reaction sets for both gas phase and 
grain surface chemistry that include extensive descriptions of 
deuterium and spin state chemistry. Our goal is not only to study 
how the abundances of light deuterated species are affected when 
grain surface chemistry is taken into account, but also to study 
spin state abundance ratios and compare these against the pre- 
dictions of the complete depletion model. 

The paper is structured as follows. In Sect.|2] we introduce 
the new chemical reaction sets and the assumptions made in con- 
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structing them. We also discuss our chemical and physical mod- 
els in general. Section [3]presents the results of our model calcu- 
lations. In Sect.|4]we discuss our results and in Sect.|5]we present 
our conclusions. 



2. Model 

We have constructed new chemical reaction sets for both gas 
phase and grain surface chemistry. The new gas phase reac- 
tion set is based on the publicly available Ohio State University 
osu_O3_2O0q3 (hereafter OSU) reaction set, which was ex- 
panded to include the spin states (i.e. ortho and para forms) of the 
light hydrogen-bearing species Hj, H2 and H^ (here, this pro- 
cess is referred to as ortho/para separation). In addition, species 
with up to 4 atoms were deuterated using a reaction cloning 
process. A new surface reaction set was also produced by per- 
forming a similar ort ho/para separa tion and deuteration of the 
surface reaction set o f ISipiial (20121). itself based on the surface 
chemistry network of lSemenov et al.l2010t the network includes 
parameters (e.g. branching ratios, activation energies if applica- 
ble) for both ionization (by either external or cosmic ray-induced 
photons) processes and reactions via the Langmuir-Hinshelwood 
mechanism (see Sect. 12.21 ). 

We present below the physical and chemical models used in 
this study and then discuss the ortho/para separation and deuter- 
ation process in detail. 



2.1. Physical model 

In this paper, we present radial abundance profiles of various 
chemical species. These profi les are produ ced identically to the 
method discussed in detail in lSipilal d2012l): we choose a modi- 
fied Bonnor-Ebert sphere (e.g. iKeto & Field! l2005t ISipila et all 
I201U ISipilal 120121) as the physical core model and combine 
chemical and radiative transfer calculations to produce self- 
consistently calculated radial profiles for both chemical abun- 
dances and temperatures (rd ust + r gas generally). We assume 
the same cooling species as in ISipilal (1201 2l) . i.e. 12 CO (and its 
isotopes I3 CO and C 18 0), C, O and O2. In this paper we concen- 
trate on comparing our new chemistry model with the complete 
depleti on model, an d for this purpose we choose the "model A" 
core of ISipilal (120121) . that is a modified Bonnor-Ebert sphere of 
mass 0.25 M Q . The density and temperature profiles of the model 
core are plotted in Fig.Q] The core has a high average density 
(«h ~ 8 x 10 5 cm -3 ) and a low average temperature {T ~ 8 K, 
depending slightly on core age), so this core should be appropri- 
ate for comparing our new results with the complete depletion 
model. The gas temperature at the core edge is different for the 
two core ages because at t — 10 5 years CO, the main coolant, 
is not yet fully depleted at the core edge (see Sioila 2012). We 
note that the gas temperature at the co re edge is slightly higher 
at t — 10 6 years than in Sipila d2012t Fig. 2). This is due to a 
slightly lower CO abundance at the core edge brought about by 
the inclusion of spin state and deuterium chemistry in the present 
model which modifies slightly t he chemistry of undeuterated 
species compared to the model of lSipilal (120121) . 

We assume that the model core is heavily shielded against 
external UV radiation - we set Ay = 10 mag at the core edge. 
The cosmic ray ionization rate is set to £ = 1.3 x 10~ 17 s -1 and 
we assume spherical grains with radius a g = 0.1 fim. The as- 
sumed initial element abundances with respect to total atomic 
hydrogen density are given in Table Q] The initial H2 ortho/para 
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Fig. 1. The density (green solid line) and temperature (black 
lines) profiles of the model core. The solid and dashed black 
lines correspond to the gas temperature at t — 10 5 and t — 10 6 
years, respectively (see also lSipilall2012l) . The dotted line repre- 
sents the dust temperature. 

Table 1. Initial elemental abundances (see text) with respect to 
total atomic hydrogen density «h- 



Species 


n(X)/n H 


H 2 (p) 


5.00 x 10-' 


H 2 (o) 


5.00 x 10~ 4 


HD 


1.60 x 10~ 5 


He 


9.00 x 10- 2 


C + 


1.20 x 10- 4 


N 


7.60 x 10- 5 


O 


2.56 x 10~ 4 


s + 


8.00 x 10~ 8 


Si + 


8.00 x 10-' 


Na + 


2.00 x 10~ 9 


Mg + 


7.00 x 10~ 9 


Fe + 


3.00 x 10- 9 


P + 


2.00 x 10" 10 


Cl + 


1.00 x 10~ 9 



ratio is (arbitrarily) set to 1.0 x 10~ 3 (we return to this issue in 
Sect. | 4 .21); the initial HD abundance is taken from Sipila et al. 
(l2010t her eafter S10) and the res t of the initial abundances are 
taken from lSemenov et all d2010l) . 



2.2. Chemical model 

The chemical code used in thi s study is essentially the same 
as the one discussed in ISipilal (1201 2l) . with some minor bug 
fixes. Gas phase and grain surface reactions are solved simul- 
taneously using the rate equation method, and gas-grain inter- 
action occurs through adsorption and thermal desorption pro- 
cesses. The rate coefficient for adsorption is A: ads = u/Scr, where 
Vj = sJdiksT gas /7mii is the thermal speed of species i, S is the 
sticking coefficient (set to unity for all species) and cr = net 
is the grain cross section. Desorption occurs mainly^ through 
transient heating of the grains by cosmic rays, with rate coef- 



See http : //www . physics . ohio-state . edu/~eric/ 



2 Thermal desorption as a separate process is also included in the 
model, but this is negligible in the temperatures considered here. 
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ficient &cr = / £td(70 K), where kju is th e thermal desorption 
rate coefficient (lHasegawa & Herbsall993h . Photodesorption is 
not included in the present model; we return briefly to this issue 
in Sect.03] 

Species adsorbed onto grain surfaces are assumed to be ph- 
ysisorbed and react via the Langmuir-Hinshelwood mechanism. 
Following the formalism of lHasegawa et al.l (Il992l) . the rate co- 
efficient for surface reactions is given by 



k ij = aK ij (R? B +Rf)/n d 



(1) 



where a is the branching ratio in the case that the reaction has 
multiple product channels, Ky is the reaction probability, Rf s 
is the diffusion rate of species i on the grain surface and rid is 
the density of grains. In addition to being destroyed in two-body 
reactions, species on the surface can be photodissociated either 
by external UV photons (A: phot = a exp(-y/Av)) or by cosmic ray 
induced UV photons (£ CR P hot = a £), although the former process 
is negligible in the high extinction environment assumed here. 

The explicit expressions of /e,y and R dlfi depend on whether 
or not quantum tunneling is assumed to occur. In case of thermal 
diffusion 



R 



diff 



diff , 



N, 



exp(-E? m /T 6mt ) 



(2) 



where v; = J2n s k B E^ / 'ifim-,; N s - n s Ana 2 g is the number of 

binding sites on the grain and E h an d Ef lS are the binding and 
diffusion energies, respectively. As in lSipiral(l2012l) . we assume 
n s - 1.5 x 10 15 cirT 2 for the surface density of binding sites and 
E dlS = 0.77 E h for the relation between diffusion and binding 
energy. If a reaction is exothermic and has no activation energy, 
Kjj — 1 , but for exothermic reactions with activation energy, /c, ; - = 
exp(-£' a /rdust) where E. d is the activation energy of the reaction. 
If quantum tunneling of light species (H and D) is assumed (see 
also Sect. l4.3t . the reaction probability and diffusion rate are - 
in the reactions involving these species - replaced by 



k u = exp[-2(a/h)(2fik B E a ) 1/2 ] 



and 



Rf« = ^exp[-2(a/h)(2mk B E?y 2 ] , 



(3) 



(4) 



respectively (lHasegawa et al.l 19921) . where /i is the reduced mass 
of the reactants. We assume a = 1 A for the barrier width. We 
note that the above assumption of a rectangular barrier is not 
valid for endothermic reactions. Some of the reactions with ac- 
tivation barriers included in the surface reaction set may be en- 
dothermic. However, because reactions with high barriers turn 
out not to influence our results significantly, we have chosen to 
omit this possib le prob l em wi th the tunneling probability. 

Contrary to ISipilal (1201 2l) . we assume 500 K for the bind- 
ing energy of H2, which is in the range of typical values as- 
sumed by other autho rs, usually somewhere betwe en 430 K 



(|Garrod&Paulvl201ll) and 600 K dCazaux et al.l2010h . In lSipilal 
(120121) . the binding energy of H2 was arbitrarily decreased to 



100K to avoid the problem of producing unphy sical amounts of 
surfac e H2 (this issue has been discussed by e.g. lGarrod & Paulvl 
l201ll) . As we now consider a binding energy value of 500 K, the 
surface H2 abundance is typically high in our models; however, 
due to the low temperature (see below) considered here, a high 
surface H2 abundance should not present a problem as the sur- 
face reactions involving H2 have high activation energies (and 
H2 is not allowed to tunnel). 



As in ISipilal (|2012[). we assume a binding energy of 1390 K 
for atomic oxygen dBergeron et al.l l2008t ICazaux et al.l 1201 lb . 
Apart from H2 and O, th e binding energie s for u ndeuterated 
species are adopted from iGarrod & Herbst (120061). Following 
the usual approach in the literature (e.g. ICazaux et al.l 120101: 
Taqu et et al.ll2013bl) . we assume that for each deuterated species, 
the binding energy is equal to that of the corresponding un- 
deuterated species. We discuss the binding energies further in 
Sect.EH 



2.3. Ortho/para separation 

To add the spin states of H}, H2 and H^ into the OSU reaction 
set, we first extracted the reactions involving these species from 
the OSU database. A Python script was written to analyze each 
reaction and to add the spin states using a predetermined set of 
rules. 

Most reactions are separated according to spin selection 
rules. As an example, consider the reaction 

Ht + O^OH + +H 2 , (5) 

with rate coefficient k. This reaction has three possible pathways: 

H+(o) + O -i» OH + + H 2 (o) 

H+(p) + oioH + +H 2 (o) (6) 

H+(p) + oXoH + + H 2 (p); 

the b ranchi ng ratios can be calculated with e.g. the method of 
Oka ([20041 ; see also Appendix [Alt. We note that in the litera- 
ture, alternative app roaches to the separation exist - for example, 
iFlower et al.l (|2006) assumed branching ratios of 1 : 2 for the sep- 
aration of HJ (p) in reactions analogous to (O, while lPagani et al.l 
(2009) assumed 1:1, as is done here. 

The above branching ratios are not applied to all reaction 
types. One example of this is charge transfer reactions where 
we have assumed that the ortho/para forms are conserved in the 
reaction. Also, whenever H2 is formed in a reaction containing 
only reactants other than HJ, H2 and H3 , such as in the reaction 



NH + + H z O — > HNO + + H 2 , 



(7) 



we assume that H2 is formed in its para state (see Appendix lAV 
This is a simplifying assumption, ensuring that we do not have 
to follow the spin states of every species with multiple protons. 
There are several special cases besides reaction (|7} - they are 
listed in Appendix lAl along with our assumptions. 

There are approximately 1000 reactions containing H^, H2 
or H^ in the OSU reaction file; performing the ortho/para sep- 
aration results in ~1700 reactions, so that the amount of new 
reactions brought about by the separation process is not signifi- 
cant (with respect to the total size of the reaction network, which 
is about 1 1600 reactions after ortho/para separation and deuter- 
ation; see Sect. l2.41 i. A similar separation process was applied to 
the surface reaction set; in this case the separation process adds 
only about 50 new reactions. 

2.4. Deuteration process 

After the ortho/para separation process was carried out, we pro- 
duced deuterated versions of both the gas phase and surface re- 
action sets. In practice, deuterons are substituted in place of pro- 
tons in each reaction and combinatorial arguments are invoked 
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Fig. 2. Radial abundances (with respect to n(H 2 )) of protons, H^ and its deuterated isotopologues and electrons (indicated in the 
figure). Panel (a) is a reproduction of the middle panel of Fig. 6 in ISipila et al.l d2010t) (showing abundances at t — 10 6 years), 
calculated with the version of the chemical code used in that paper (see text). The two middle panels use the same physical model 
as in panel (a), but chemical evolution has been calculated with the present version of the chemical code (see text) with quantum 
tunneling turned either on (panel b) or off (panel c). Panel (d) plots again the same physical situation but calculated with the new 
reaction network presented in this paper, with initial heavy element abundances set to zero. 



(as in the case of ortho/para separation) to work out branching 
ratios, where applicable. For example, deuterating the reaction 



NH + + H 2 — > NH+ + OH 



(8) 



results in 9 new reactions including, for example, the reactions 



NH + + HDO — -> NHD + + OH 



NH + + HDO- 



H 



• NH+ + OD , 



(9) 



where complete scrambling, i.e. a statistical approach, is as- 
sumed. Our de uteration routine is similar to that presented in 
iRodgers & Millar! (1 1996b - recently an analo gous method has 
also b een applied to the os u_2QQ9 network by lAlbertsson et all 
(120111) . As pointed out by iRodgers & Millar! (119961) . the com- 
plete scrambling assumption is not appropriate for all reactions, 
particularly for those involving more complex molecules with 
OH endgroups, for example. However, as we are interested in 
the chemistry of relatively simple species in the present paper, 
we have chosen a rather relaxed approach toward this problem 
and assume that the complete scrambling assumption is gener- 
ally valid. Of course, there are a number of special cases; for 
example, in charge transfer reactions no atoms are interchanged, 
and the deuteration routine is designed to make sure this does not 
occur in the deuterated version of the reaction either. Other im- 
portant special cases are deuterated analogues of reactions such 
as 



HJ + H 2 



H?0 + + H 2 



(10) 



We have assumed that reactions such as this proceed by H^ do- 
nating a proton and thus, for example, the reaction D 2 H + + H 2 
cannot result in HD 2 + + H 2 . 

Deuterated reactions involving only the lightest species (up 
to helium) are taken from previous works. The deuteration chem- 
istry of the light species in the context of a complete depletion 
model (no heavy elements in the gas phase, and no surface chem- 
istry except f or the formation o f H 2 , HP a nd D?) was discusse d 
originally by IWalmslev et all (|2004|) and iFlower et all (|2004|) . 
and subsequently by many authors. The goal of the present paper 



is to study how the inclusion of heavy elements and a descrip- 
tion of gas-grain interaction including surface chemistry affects 
the gas-phase chemistry of H3 and its deuterated isotopologues; 
we have therefore incorporated the reaction set used in S10 into 
the expanded OSU set so that the final reaction set is consistent 
with the complete depletion model as far as the chemistry of 
light species is concerned. In practice, we replaced the appropri- 
ate reactions in the new reaction set with their S10 counterparts 
and added any reactions from S10 that did not exist in the new 
modejj. This allows for a maximally consistent comparison be- 
tween the two cases. 

We have similarly added to the final (deuterated and spin 
state-separated) surface reaction set the reactions a nd the associ- 
ated a ctiva tion energies include d in the models of ICazaux et al.l 
(120101) and lCazaux etail d201ll) . 

As a final step, all reactions in the final deuterated gas phase 
and surface reactions sets were checked for elemental balance, 
i.e. that both sides of all reactions contain the same amount of 
elements. The final gas phase reaction set, including deuterium 
and spin state chemistry, contains about 11600 reactions. The 
final surface reaction set contains about 1350 reactions. 



3. Results 

3. 1. Comparison of the new chemical code against the S10 
code 

A previous (gas-phase chemistry only) version of our chemical 
code was used in S 10 to study deuterium chemistry in the com- 
plete depletion limit; to test the new reaction set discussed here, 
and the chemical code itself, we reproduced some of the results 
of S10 using both the old and the new versions of the chemical 
code (Fig.[2j. Panel (a) is a reproduction of the middle panel in 
Fig. 6 in S10 (abundances are shown at t — 10 6 years), calcu- 
lated using the old version of the chemical code and adopting 
the same density and temperature profiles, the same reaction set 



3 For example, the OSU set does not include the reaction HJ + H 2 — > 
HJ + H2, but when considering spin states (and deuterated forms) ex- 
plicitly, the variants of this reaction need to be a dded to the model. This 
important system has been recently analyzed bv lHugo et al.l(l2009l) . 
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and assuming the same values for the physical parameters as in 
S10. It should be noted that because only gas-phase chemistry 
was considered in S10, the grain-surface formation mechanisms 
for H2, HD and D2 were hard-coded into the program and no sur- 
face species as such were included. The new version of the code, 
however, inc ludes surface chemistry as described above and in 
ISipilal (120121) : to compare these two approaches, we plot in pan- 
els (b) and (c) the results of calculations otherwise identical to 
those of panel (a), but calculated using the new code with quan- 
tum tunneling either included (panel b) or excluded (panel c). It 
can be seen that without quantum tunneling, the results of S10 
are not reproduced. The reason for this is the very low tempera- 
ture (~6.5 K in the center) of the model core; without quantum 
tunneling, surface production of HD is inefficient at this tem- 
perature and eventually leads to a drop in gas phase HD abun- 
dance. This in turn increases the H3 and H2D + abundances and 
decreases the D2H + abundance, while the D^ abundance is rel- 
ative ly unaltered. The re a sons for these effects are discussed in, 
e.g.. IFlower etal.l (120041) . iFlower et all (120061) and IPagani et al] 
d2009t) . With quantum tunneling included, the new code gives 
virtually identical results compared to the old code because of 
the rapid conversion of surface H and D to HD. It should be 
noted that the creation rate of HD on gr ain surfaces depends crit - 
ically on the assumed temperature (e.g. lCazaux & Spaansl2009l) ; 
already at ~8 K, the difference between the tunneling and no tun- 
neling cases is very small as demonstrated by comparing results 
at radii > 2000 AU in panels (b) and (c) of Fig.[2l where the 
temperature is ~ 7.5 K. 

In panel (d) we plot again the same physical situation as in 
panels (a-c), using the new reaction set described above but with 
initial heavy element abundances set to zero (and quantum tun- 
neling included). As would be expected since the light element 
chemistry originates mainly from the S10 set (see Sect. 12.4) . the 
abundances are practically identical to those presented in pan- 
els (a) and (b). Even though the deuteration process described 
above introduces some reactions involving only the light species 
that do not exist in the S10 set, it turns out that these reactions 
(such as some reactions involving the H and D anions) are 
rather insignificant (at least for this set of physical conditions). 

In what follows, all modeling results correspond to the new 
ortho/para separated and deuterated reaction set with quantum 
tunneling switched on, unless otherwise stated. 

3.2. Comparison with the complete depletion model 

In the complete depletion model, it is assumed that "heavy" 
species (i.e. those containing elements heavier than He) are 
frozen onto grain surfaces at high gas densities, allowing the 
deuteration chemistry of H3 to proceed unhindered. The deuter- 
ation chairi_cTJTJj^higWy dependent on the abundance of HD 
(e.g. IWalmslev et al.l 12004]) . which is produced on grain sur- 
faces and is thus dependent on the available amount of atomic 
H and D. In the complete depletion model, grain-surface atomic 
H and D are spent only in the reactions creating H2, HD or D2; 
in essence, it is assumed that heavy species on grain surfaces 
are locked in unreactive species, such as water or methanol. 
However, in a model where heavy elements are not totally de- 
pleted, the large amount of reacting species on the surface may 
reduce the amount of HD produced because of the additional re- 
action pathways for H and D (compared to the complete deple- 
tion model), and thus affect the abundances of HJ and its deuter- 
ated isotopologues in the gas phase. 

To demonstrate this, we plot in Fig.|3]the radial abundances 
of Ht and its deuterated isotopologs using the reaction set of 



either S10 (panel a) or this work (panels b and c); also plotted 
are the protons and electrons. Panels (a) and (c) correspond to 
a core age of 10 6 years and panel (b) corresponds to t - 10 5 
years. It can be immediately seen that the abundances change 
dramatically depending on the adopted reaction set and that in 
the present model the abundances are highly time-dependent 
even at advanced core ages (in the complete depletion model, the 
abundances reach steady-state across the core at about t — 2x 10 5 
years). In the dense center of the core, where one would expect 
a high degree of d euteration of Ht ba s ed on the complet e de- 
pletion model (e.g. IFlower et all 12004 IPagani et ai1l2009l) . the 
present model predicts a completely opposite situation at t — 10 6 
years - in this case, the degree of deuteration increases toward 
the edge of the core, i.e. toward lower density. This can be under- 
stood through the surface chemistry; we plot in Fig.|4]the abun- 
dances of HD, HDO* (in this paper, an asterisk represents a sur- 
face species) and H^ and its deuterated isotopologs as a func- 
tion of time in the innermost (left panel) and outermost (right 
panel) shells of the model core. In the dense center of the core 
(«h ~ 10 6 cm 3 ), H* and D* react preferentially with O* to pro- 
duce OH* and OD*, which react again with H* and D* to produce 
H2O* and HDO* (these processes are efficient because oxygen is 
initially in atomic form; see Tableland Sect. 14. U . This process 
almost completely suppresses the formation of HD* through the 
reaction of H* with D* (see also Sect. l4.31 >. Thus, since its abun- 
dance is not efficiently regenerated on grains, HD ultimately de- 
pletes from the gas phase. HD depletion coincides with heavy 
element depletion - this is because the H3 abundance increases 
with heavy element depletion and starts to destroy HD, which 
is not efficiently created on grain surfaces owing to the reasons 
mentioned above. 

The HD abundance experiences a slight increase at about t = 
3.5xl0 5 years. This arises because OH* and OD* disappear from 
the grain surfaces as they are processed into H2O* and HDO*; 
this frees up some H* to react with HDS* and HDCO* to produce 
HS* + HD* and HCO* + HD*, respectively (see also Sect.l43T>. 
The HD* thus formed desorbs and momentarily regenerates the 
HD abundance. 

At the lower density of the core edge (~ 10 5 cm -3 ), repre- 
sented by the right panel of Fig.|4] HD depletes less than at very 
high density (~ 10 6 cm' 3 ), implying a higher deuteration degree 
of H^ at an advanced core age. The slower depletion of HD at 
the lower density is due to the overall depletion timescale being 
longer in these conditions, which suppresses deuterium chem- 
istry (as seen in the behaviour of H^ deuteration at early times), 
and thus there is somewhat less D* available to form OD* (and 
HDO* through OH* + D*). As a consequence, the timescale for 
deuterium incorporation into HDO* is longer than at high den- 
sity. Also in this case, total HD depletion is prevented by the 
reactions H* + HDS* and H* + HDCO*. 

Figures [3] and [4] demonstrate that the deuteration degree of 
H^ is highly dependent on time. Indeed, at t — 10 5 years, the 
agreement of the present model with the complete depletion 
model is much better than at t = 10 6 years. With the exception 
of DJ , for which we obtain much lower abundances than in the 
complete depletion model, the agreement of the two models is 
rather good in the core center at t = 10 5 years, but less so at the 
core edge where depletion is not as efficient. The detailed behav- 
ior of the deuteration degree depends on many factors, including 
the binding energies of various species, but the main result of 
Figs.[3]and|4]is the depletion of HD which does not occur in the 
complete depletion model. 

From an observational standpoint it is interesting that the ra- 
dial abundances of H2D + and D2H + are similar to each other in 
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Fig. 3. Radial abundances (with respect to «h) of protons, electrons and H^ and its deuterated isotopologs (indicated in the figure). 
Panel (a) corresponds to the complete depletion model at t — 10 6 years (using the reaction set of S 10), panels (b) and (c) correspond 
to the present model at t — 10 5 years and t — 10 6 years, respectively. 




Fig. 4. Abundances (with respect to «h) of selected species (indicated in the figure) as functions of time in the innermost shell 
(R ~ 240 AU, « H ~ 1-9 x 10 6 cirT 3 ; left panel) and the outermost shell (R ~ 3600 AU, n H ~ 2.2 x 10 5 cirT 3 ; right panel) of the 
model core. 
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Fig. 5. Nuclear spin state abundance ratios of H3 and its deuterated isotopologs and of JU (multiplied by 10 5 ). Panel (a) corresponds 
to the complete depletion model at t — 10 6 years, panel (b) to the present model at t — 10 5 years and panel (c) to the present model 
at t — 10 6 years. 
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the new model at t — 10 5 years (see Fig. [3), although their abun- 
dances are clearly lower in the new model than in the complete 
depletion model at the core edge. This agrees with observational 
results that the two species are present with comparable abun- 
dances dParise et alj|201 UlVastel et alj|2012l) . 

It is also observed in Fig. [3] that the ionization degree de- 
creases when one adds heavy elements to the model. This is be- 
cause in the present model, electrons are mainly removed in dis- 
sociative recombination reactions with HCO + and HCNH + ; the 
rate coefficients of these reactions are about a factor of 3 higher 
than those of H^ and its deuterated forms, leading to a lower ion- 
ization degree in a model with heavy elements when compared 
to the one with zero heavy element abundances. 

3.3. Ortho/para ratios 

In Figs. |2] E] and |U we have plotted the abundances of HJ and 
its deuterated isotopologs as sums of the abundances of their 
respective nuclear spin states. To illustrate how the ortho/para 
ratios of these species differ in the present model compared with 
the complete depletion model, we plot in Fig. [5] the ortho/para 
ratios of H+, H 2 D + , D 2 H + and H 2 (multiplied by 10 5 ) and the 
meta/ortho ratio of D^ as predicted by both the complete deple- 
tion model and the new model presented here. The panels de- 
pict the same core ages as in Fig. [3] Figure [5] is supplemented 
by Fig. [6] where the spin state ratios are plotted as a function of 
time in the innermost and outermost shells of the model core. 
As could be expected based on the evolution of the HD abun- 
dance in the present model, the spin state abundance ratios of 
the variou s species, which depend on the abundances of HD and 
ortho -H? (IWalmslev et al.ll2004l iFlower et al.ll2004t ISipila et all 
1201 Ol) . are somewhat different in the two models, depending also 
on core age. The H^ , H 2 D + and D 2 H + spin state ratios are sim- 
ilar (within a factor of a few) in both models, particularly at 
t — 10 5 years, but there are larger differences in the behavior of 
the DJ (m/o) ratio and the H 2 (o/p) ratio depending on the model. 

One can see from Figs. [5] and [6] that the o/p ratios (and the 
spin temperatures, r sp j n ) of H 2 , H^ , and H 2 D + decrease as the 
core evolves. The effect is particularly marked for H 2 , with r sp i n 
dropping below 15 K after 10 6 years in both the innermost and 
the outermost shells of the model core. Also for D 2 H + and D3 , a 
slight decrease in r sp ; n (i.e. an increase in the o/p and m/o ratios, 
respectively) can be seen at late times (after 10 6 years) when 
compared with early times (< 10 4 years). The D^m/o) ratio 
settles at ~ 1 which indicates a clearly super-thermal r sp i n of 
42 K, but the ratio experiences a peak at around 10 5 years with 
a m/o ratio of ~ 8 (r spin ~ 16 K). For H+, H 2 D + and D 2 H + , 
the spin temperatures settle between 15 and 30 K at late stages 
of core evolution. 

The D} (m/o) ratio is much lower in the core center at t — 10 6 
years in the present model than in the complete depletion model. 
At the high density of the core center, the ratio first begins to in- 
crease following heavy element depletion (initiating the deuter- 
ation chemistry), but then decreases again due to HD depletion - 
comparing Figs. [4] and |6]reveals that a decrease of the D^ (m/o) 
ratio indeed coincides with HD depletion. This is because meta- 
Do is most efficiently created in a spin-state-conversion reaction 
of ortho-DJ with HD and in the reaction between D 2 H + (o) and 
HD; therefore HD depletion leads to a decrease of the D^ (m/o) 
ratio. 

The H 2 (o/p) ratio is very different in the present model than 
in the complete depletion model. The effect of HD depletion can 



be seen here as well; H 2 (o) is most efficiently created in the gas 
phase in the reactions 



H+ + HD — » H 2 D + + H 2 , 
H 2 D + + HD — > D 2 H + + H 2 , 



and 



D 2 H + + HD — > D+ + H 2 



(ID 
(12) 

(13) 



so that HD depletion directly affects H 2 (o) creation. HD deple- 
tion is stronger and occurs earlier at the core center than at the 
edge, and this accounts for the differences seen between panels 
(b) and (c) in Fig. [5] Note that the gas phase production of H 2 (o) 
is now the most important factor controlling the H 2 (o/p) ratio 
because quantum tunneling is included; see Sect. 14.31 

Interestingly, the present model predicts very similar H 2 D + 
and D 2 H + (o/p) ratios at t — 10 5 years compared with the com- 
plete depletion model even at the core edge, where the H 2 D + and 
D 2 H + (ortho+para) abundances differ from the complete deple- 
tion model (Fig.[3}. Finally we note that the ortho/para ratios 
predicted by our models at t — 10 6 years are for the mo st part 
compatible with recent modeling results of Vastel et al. (1201 2t 
their complete depletion case), although our chemical network 
is much larger and the present model includes a time-dependent 
treatment of depletion. 

4. Discussion 

4. 1 . Initial heavy element abundances 

In this paper, we assume that the gas is initially atomic, with the 
exceptions of hydrogen being in molecular form and deuterium 
being locked in HD (see Table[T]i. In other words, we do not take 
into account any possible chemical history of the core, i.e. the 
possibility that some of the heavier elements could be locked in 
ices at the start of the calculation. This can be particularly impor- 
tant for HDO* which is created from OH* and OD* and is thus 
dependent on the amount of available atomic oxygen. The as- 
sumed initial HD abundance corresponds roughly to the average 
gas phase D/H ab undance ratio, (D/H) gas , in the Local Bubble 
(I Wood et ai1l2004l) . Substantially lower values of (D/H) gas have 
been determi ned towards lines of sight with large column den- 
sities. iLinskv et al.l d2006l) find a correlation between deuterium 
depletion and the depletion of metals, and suggest that D can be 
incorporated in carbonaceous dust grains. The assumption about 
the initial gas phase deuterium abundance determines the over- 
all degree of deuteration, but is not likely to affect the relative 
abundances of deuterated species. 

Because the present model does not include a description 
of core collapse (from, e.g., a diffuse cloud into the dense 
configuration that is the present model core), we have investi- 
gated the effect of the initial abundances on HD depletion by 
first calculating chemical evolution in a low-density cloud with 
«h = 2 x 10 3 cm~ 3 , visual extinction Ay - 3 mag and tem- 
perature r gas = rd ust = 10 K (which should be appropriate for 
these values of density and visual extinction, luve la & Ysardl 
201 1). The abundances of all species corresponding to three dif- 
ferent times (10 4 , 10 5 and 10 6 years) were then extracted and 
used as initial abundances for the dense core. Figure Q shows 
the result of such calculations in a single-point chemical model 
with «h = 2 x 10 6 crrT 3 , temperature r gas = Td ust - 7.5 K 
and Ay - 100 mag, corresponding roughly to the interior of 
the model core discussed in this paper. It is observed that the 
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Fig. 7. The abundances of HD and HDO* as a function of time 
at a density of «h = 2 x 10 6 cirT 3 calculated for progenitor cloud 
(see text) ages of 10 4 (solid lines), 10 5 (dashed lines) or 10 6 (dot- 
ted lines) years. 



older the progenitor cloud, the longer it takes for deuterium to 
get locked into HDO*. This is because after 10 6 years of progen- 
itor cloud evolution, about a half of the initial oxygen abundance 
has been converted to H2O*, which limits somewhat the surface 
reactions leading to the production of HDO* in the dense core 
phase (the HDO7H2O* ratio is only about 10 3 at the end of 
the progenitor cloud evolution). Notably, very little HD deple- 
tion takes place in the diffuse cloud. We also tested the effect 
of more extreme initial abundances by locking all oxygen ini- 
tially in H z O* (50%) and CO* (50%) and all nitrogen in NH; - 
in these calculations, HD still depletes (by about two orders of 
magnitude in ~ 10 6 years) because surface D gets locked into 
multiple singly deuterated surface species (NH2D*, HDCO* and 
HDO*). 

While the details of the chemistry change when one starts 
from molecular initial abundances (with respect to the model 



with initially atomic heavy element abundances), our main con- 
clusion (HD depletion) remains unaffected. Figure [7] demon- 
strates this for «h = 2 x 10 6 crrT 3 , but we have also performed 
test calculations at lower densities to confirm that HD depletion 
is largely unaffected in those conditions as well. 



4.2. Initial H 2 ortho/para ratio 

Dense cores condense from diffuse molecular clouds where the 
observationally derived ortho/para column density ratios of H 2 
indicate s pin temperatures in th e range 50-70 K (H 2 (o/p) ~ 
0.3 - 0.8, ICrabtreeetaDl201ll) . ICrabtree et aD (1201 lb showed 
that the H 2 spin temperatures in diffuse clouds correspond to 
the gas kinetic temperatures because the ortho/para ratio of H 2 
is thermalized through collisions with H + . Also in dense clouds 
with elevated temperatures, the spin-state conversion reactions 



H + + H 2 (o) -» H + + H 2 (p) 

and 

H + + H 2 (p) -» H + + H 2 (o) 



(14) 



(15) 



govern the H 2 (o/p) ratio. As discussed by Flower et al. {2006; 
their Sect. 3.1), the spin temperature of H 2 follows closely the 
kinetic temperature down to T^ n ~ 20 K. The H 2 (o/p) ratio cor- 
responding to this temperature is 1.8 x 10 -3 . Below 20 K, the 
endothermic para-ortho conversion reaction with H + (reaction 
IT~5T > starts to lose importance, and H 2 (o) production is dominated 
by formation on grains. The consequence is that with a decreas- 
ing kinetic temperature, the H 2 spin temperature rises above the 
kinetic temperature. In cold, dense cores {T^ n ~ 10 K), the com- 
petition between the more favored ortho production on grains 
and the efficient para-ortho conversion in reactions with H + and 
H| (and its deuterated counterparts) in the gas phase leads to an 
H 2 spin temperature which lies a few K above the kinetic temper- 
ature. This point is illustrated in Fig.|8]at high values of gas den- 
sity («(H 2 ) = 10 5 crrT 3 ) and visual extinction (Ay = 10 mag); 
the gas (solid line) and dust (dashed line) H 2 spin temperatures 
lie clearly above the kinetic temperature (dotted line) at low val- 
ues of 7\j n . 
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Fig. 8. The spin temperatures of gas-phase (solid line) and grain- 
surface (dashed line) H2 as a function of the kinetic temperature 
forn(H 2 ) = 10 5 cnT 3 andA v = 10 mag. 



Based on the ideas presented above we have chosen, some- 
what arbitrarily, the initial value 1 .0 x 1CT 3 for the H2(o/p) ratio 
in our chemistry model (see Table [TJ. This ratio corresponds to 
an H2 spin temperature of 18.7 K. Although the inital ortho/para 
ratio depends on the thermal history of the core, the spin tem- 
perature should in any case lie somewhere between the current 
kinetic temperature and 20 K, because during core condensation 
and cooling, H2 has probably been thermalized down to 20 K. 
Therefore, we think the initial value quoted above is more real- 
istic than the thermal value at 10 K (3.6 x 10~ 7 ) or the statistical 
ratio at high temperatures (3) which is the assumed ortho/para 
formation ratio on grains. The choice of the initial abundances 
has some effect during the first ~ 10 5 years in our simulations; 
at late times they are forgotten because of the high density of the 
model core. 



4.3. Quantum tunneling 

Whether diffusion of species via quantum tunneling occurs on 
grain surfaces or not is still a subject of some debate in modern 
astrochemistry. There is so me experimental evidence that tun- 
neling does not occur (e.g. iKatz et aD 1 19991) and on the other 
hand some modeling r esults indicate that tunne ling should be 
taken into account (e.g. lCazaux & Tielensll2004l) . In this paper, 
we have included quantum tunneling so that our analysis would 
be maximally consistent with the complete depletion model (see 
also Sect. l3.1l and Fig.|2j». However, it is also prudent to take into 
account the possibility that tunneling does not occur and to see 
whether HD depletion is affected by this assumption; with this in 
mind, we have carried out calculations similar to those presented 
in the earlier sections with quantum tunneling turned off. The re- 
sults of these calculations are presented in Fig. [9] which should 
be compared with Fig. [4] It is observed that while the detailed 
behavior of the abundances as a function of time changes when 
quantum tunneling is turned off, the main result of this work - 
the ultimate depletion of HD at high density - is not affected. 

The details of the (deuterium) chemistry change when quan- 
tum tunneling is turned off, which can be explained as follows. 
When one assumes quantum tunneling, many of the surface reac- 
tion pathways for H and D that would be otherwise unavailable 
due to high activation energy become accessible because tun- 



neling increases the reaction probability and the effective rate 
coefficient by increasing the diffusion rate (see Sect. 12. 2b . This 
effectively reduces the abundances of H* and D* so that it be- 
comes increasingly difficult to form H*. and HD* through the 
basic reactions H* + H* and H* + D*, respectively. Indeed, 
the main pathway for producing HD* at advanced core ages is 
H* + HDCO* (as discussed in Sect. l3.2b - when tunneling is 
included. However, without tunneling, all of the channels for 
which E a + are effectively closed off owing to the very low 
temperature, and HD* is mainly produced in H* + D* after OH* 
and OD* are processed into H2O* and HDOI3 This reaction is 
more efficient than H* + HDCO* in the tunneling case and thus 
gives a slightly higher HD* yield, and for this reason HD* de- 
pletes slightly slower in the model with no tunneling (also ob- 
served in the form of a somewhat more pronounced bump in the 
HD abundance at high density); the decreased HD depletion also 
slightly increases the abundances of the deuterated isotopologs 
of H3 at late times. 

The model without tunneling also yields a larger amount of 
H2(o) at late times compared to the model with tunneling, which 
is demonstrated in Fig.[10](to be compared with Fig. [6}. The rea- 
son for this is again the increased H* abundance which leads 
to efficient H2(o) creation through the H* + H* reaction. When 
tunneling is included, H2(o) is not efficiently created on grain 
surfaces because 1) the H* + H* reaction is suppressed, and 2) 
other surface reactions preferentially create H2(p) owing to the 
adopted ortho/para separation rules (see Appendix |A}. An ob- 
servationally important consequence is the increased H2D + (o/p) 
ratio at late times when tunneling is turned off. 

To summarize, the inclusion or exclusion of quantum tun- 
neling modifies the detailed behavior of the deuterium chemistry 
and particularly the spin state ratios, but the main result of this 
paper (eventual HD depletion at high density) is not affected. 
While our assumptions regarding the production of H2(o) natu- 
rally influence its abundance in the model, more (experimental) 
evidence is needed to justify either including or excluding quan- 
tum tunneling on grain surfaces so that the detailed behavior of 
the spin state abundances at advanced core ages can be better 
constrained. 



4.4. Binding energies 

Experimental evidence suggests that the binding energies of 
light species on gr ain surfaces depend st r ongly on the proper - 
ties of the ice (e.g. iHornekaer et al.ll2005t lAmiaudet al.[|2006l) . 
Because in a realistic situation the properties of the surface 
(porosity, content of the ice) change as a function of time, one 
should consider time-dependent binding energies for the various 
species instead of the constant values considered here (and in 
the majority of other works). Without an appropriate model, it 
is difficult to quantify whether HD depletion would be affected 
if this were done, but to study the influence of varying binding 
energies we have run some test calculations using a single-point 
chemical model corresponding to the conditions in the center of 
the model core, varying the binding energies of H and H2 and 
their deuterated forms between 300 K and 600 K (which should 
cover the typical values on either non -porous or porous amor- 
phous solid water; iTaquet et al.ll2013bl) . In all test calculations, 
HD eventually depletes from the gas phase, however the deple- 



4 Since the H* + O* and D* + O* reactions are barrierless, the inclusion 
of tunneling does not change the reaction probability, but does increase 
the effective rate coefficient. However, these reactions are efficient both 
with and without tunneling owing to the large O* abundance. 
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tion is somewhat less pronounced with increasing H and D bind- 
ing energy when tunneling is not included (because a higher H 
or D binding energy translates to lower surface reaction rates so 
that HDO* is not formed as efficiently). If tunneling is included, 
the variation in binding energy does not influence HD depletion 
in the range of binding energies studied (unless H2 is allowed to 
tunnel, see Sect. 14.: 



In this work, we assume that the binding energy of each 
deuterated species is the same as that of the undeuterated ana- 
logue. This is an approximation since in principle, one would 
expect deuterated species to be somewhat more strongly bound 
to grain surfaces than their undeuterated c ounterparts due to 
their slightly higher mass dTielensI 119831: iPerets et all l2005t 
iKristensen et al.ll201 ll) . However, we do not expect taking the 
binding energy difference between non-deuterated and deuter- 
ated species into account to influence the main result (HD de- 
pletion) of this paper; if we assumed higher binding energies for 
deuterated species, this would probably only lead to more effi- 
cient HD depletion owing to decreased desorption efficiency. 



4.5. The HDO* abundance 

Looking at Fig. [4] it seems that the conversion of HD into HDO* 
depends only very loosely on density. In these high-density 
cases, oxygen is finally locked mainly in H2O* while deuterium 
is locked in HDO*. However, at lower densities, the recycling 
of deuterium is slower and the late-time abundance of HDO* 
is consequently lower than at high density. This is illustrated in 
the left panel of Fig.QT| where the ratio of HDO* to the ini- 
tial gas phase HD abundance at different densities is plotted 
as a function of time. For these calculations, we have assumed 
T g -ds - Tdust = 10 K and Ay = 10 mag. Evidently, the relative 
amount of deuterium locked in HDO* at low density is small 
even at advanced core ages. Because HDO* is the most abun- 
dant D-bearing species on grain surfaces in our model, it follows 
that HD is hardly depleted at all at low density. 

In the right panel of Fig.[TTJ we plot the HDO*/H 2 0* ra- 
tio at different densities as a function of time. This plot can be 
compared with other studies of water ice deuteration, such as the 
recent study by iTaquet et al.l d2013bl) : although our model is in 
many respects different from theirs, we predict similar values for 
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Fig. 11. Left: The ratio of the HDO* abundance to the initial gas phase HD abundance at different densities (indicated in the figure) 
as a function of time. Right: The HDOVH2O* ratio at different densities as a function of time. 



the HDO7H2O* rati o for typical core lif etimes (10 5 - 10 6 years; 
see Figs. 8 and 10 in lTaquet et al.l2013bl) . The agreement is very 
good consid ering the diff erences between our model and theirs; 
for instance jTaquet et al.l d2013bl) consider fixed values of the H2 
ortho/para ratio, whereas in our model the ratio is calculated as a 
function of time. Furthermore, in our model the binding energies 
of the various species do not vary with the H2 content on the sur- 
face (see Sect. 12.21 ), which should make some difference on the 
progression of the surface chemistry, although this is difficult to 
quantify without a direct comparison of the two models (see also 
Sect.l4Tb. Finally, our model (unlike that of iTaquet et~ail20 1 3bl) 
does not include a treatment of photodesorption. Consequently, 
the desorption efficiency is somewhat lower in our model which 
may translate, in particular, to slightly hig her HDO* abundance s 
in the present model than in the model of ITaquet et al.l (l2013bl) . 
It should be noted that it is uncle ar if HD dep l etion oc curs in 
a multilayer model, such as that of ITaquet et al.l d2013bl) . Based 
on our results, HD depletion occurs because deuterium is mostly 
locked in HDO*, i.e. the late-time HDO* a bundance is similar t o 
the initial HD abundance. In the model of lTaquet et al.l d2013bl) . 
the H2O* and HDO* abundances are about an order of magnitude 
lower than in our model in similar physical conditions, so that 
the entire deuterium reservoir in HD cannot be transferred solely 
to HDO*. However, in a multilayer model deuterium might still 
be efficiently locked into multiple surface species (e.g. H2O, 
H2CO, NH3, ...) which could influence the gas-phase HD abun- 
dance - this issue should be investigated. 



lAikawa et al.l d2012l) have recently studied deuterium chem- 
istry in the context of a collapse model, and they predict 
HDO*/H 2 0* ratios of ~ 0.01. Their model does not in- 
clude spin state chemistry, and thus gives an upper limit on 
the str ength of deuterium fractionation (whic h is reduced by 
H>(o1:lFlower et alJl2006t ITaquet et al.ll2013bl) . In the model of 



lAikawa et al.l d2012l) . the sticking coefficient is set to 0.5 whereas 
we use a value of unity, which means that the relative amount of 
reactive species on the surface should be l arger in our model a t 
any given time. Finally, in the model of lAikawa et al.l (120121) . 
the model core spends a relatively short time in the prestellar 
phase (the protostar is born in 2.5 x 10 5 years) - comparing the 



HDOVH2O* ratio given by their model near the co re edge where 
wh ~ 10 5 cm" 3 and T ~ 10 K (Figs. 1 and 3 in lAikawa et alJ 
l2012h with our results (Fig.QT), the agreement is again rather 
good. 



ICazaux et al.l d2011l) reported modeling results yielding very 
low HDOVH2O* ratios at low temperatures (and a strong tem- 
peratu re dependence of the ratio). In the model of ICazaux et al.l 
d2011l) . H2 is allowed to tunnel and this allows H2O* to be ef- 
ficiently produced in the O* + H*, reaction des pite its high ac- 
tivation barrier (3000 K in ICazaux et al.l 1201 ll) . which in their 
model leads to a low HDOVH2O* ratio at low temperatures. 
We have performed test calculations including the reactions of 
ICazaux et al.l d2011l) in the chemical model and allowing H2 to 
tunnel as well (in addition to H and D) - the results of these cal- 
culations are plotted in Fig. [12] which assumes the same physi- 
cal parameters as in Fig.QT] Evidently, HDO* is produced less 
efficiently at increasing densities. This is because H2 tunnel- 
ing allows the barriers associated with the reactions involving 
H2 to be more easily overcome, and the O* + H*, reaction (in 
particular) becomes important. Consequently, O* and OH* are 
used up fast and the HDOVH2O* ratio levels off before the gas- 
phase deuterium chemistry, releasing atomic D, takes over. Thus, 
we find low HDOVH2O* ratios (< 10 3 ) even at high density, 
which leads to little HD depletion. We note that while this is 
an interesting result, more experimental data not only on tun- 
neling on grain surfaces but also on the O* + H*, r eaction at 
low te mperatures is needed to investigate this issue. lOba et all 
d2012l) claim to find no evidence for this reaction to proceed 
at low temperature. However, because of the difficulty to carry 
out laboratory experiments with atomic oxygen, more experi- 
ments are needed to confirm the Oba et al. result. Finally, the 
results of Fig. [12] are also affected by the H2 binding energy; 
for example for £^ = 300 K, surface H2 desorbs fast which 
greatly decreases the rates of the reactions involving H2, and we 
get results closer to Fig.QT] Therefore, a model including time- 
dependent binding energies, possibly yielding less H*, than the 
static-binding-energy model considered here, could also yield 
high HDO* abundances at high density even when H2 tunnel- 
ing is included. This point warrants further investigation. 
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Fig. 12. As Fig.Qj] but assuming H2 tunneling (in addition to H and D). 



4.6. Observational constraints on the HDO* abundance 

The depletion of heavy elements and the accompanying increase 
in the atomic D/H abundance ratio in the gas phase is expected 
to result in high degrees of deuteration in water, formal dehyde, 
and m ethanol incorporated in the icy mantles of grains (iTielensI 
Il983l) . Observations suggest, however, that deuteration is clearly 
more pronounced in H2CO* and CH3OH* than in H2O* (see e.g. 
iTaquet et al.ll2013bl and the references therein). This has been 
explained by the fact that water ice is formed at an early stage 
when CO depletion is not marked. Observations of the O-H and 
O-D streching bands in the infrared have provided upper limits 
of ~ 1% for the solid HDO*/ H?0* ratio towards intermediate- 
mass and low-mass protostars (iParise et al.ll2003t Partois et al 
120031). The O-D band at 4.1 yum is weak and broad. iGalvez et al 
d2011l) estimate that the detection of HDO in amorphous ice is 
nearly impossible if the HDO7H2O* abundance ratio is less than 
a few percent. 

The gas-phase HDO/H2O ratios in hot corinos around 
solar type protostars, determined recently using interfero- 
metric observations, range from ~ 1 0~ 3 to a few per- 
cent (|Jjz>rg ensen & van Dishoeck 2010t iPersson et al.l 1201 3t 
IVisser et aLll20 13: Taqu et et al1l2013al) . These ratios, which are 
likely to reflect the deuteration of water ice before desorption, 
conform with relatively low HDOVH2O* ratios inferred from 
infrared absorption measurements. 

The HDOVH2O* column density ratio towards the center of 
our model core is ~ 0.03 and ~ 0.07 at t — 10 5 yr and t - 10 6 yr, 
respectively. These values are close to the high end of the range 
of gas-phase HDO/H2O abundance ratios derived in hot cori- 
nos. In principle, the HDO* column density in the model core 
becomes large enough to be determine d through infra r ed abs orp- 
tion. Adopting the parameters used by iDartois et al.l (120031) . we 
estimate that the peak optical thickness of the 4. 1 /mi O-D ab- 
sorption band from amorphous HDO ice (FWHM ~ 0.2 //m) is 
about 0.02 and 0.18 at the times t - 10 5 yr and t - 10 6 yr, re- 
spectively. These values are given for an offset from the centre 
corresponding to one half of the outer radius where the visual 
extinction through the core is Ay ~ 23 mag. The HDOVH2O* 
abundance ratios at the two times quoted are ~ 0.02 and ~ 0.06. 
The inclusion of a low-density envelope decreases these values 



and increases the extinction through the cloud, but the determi- 
nation of the ratio may be feasible towards a background star in a 
situation corresponding to a late evolutionary stage of the model 
core. 



5. Conclusions 

We have studied deuterium chemistry in starless cores with a 
gas-grain chemical model utilizing new chemical reaction sets 
for both gas phase and grain surface chemistry that include the 
spin states of the light species (Hi, H2 and Ht ) and deuterated 
forms of species with up to 4 atoms. It was found that at the high 
densities («h ^ 10 5 cm" 3 ) and low temperatures (< 10 K) ap- 
propriate to the centers of starless cores, HD eventually depletes 
from the gas phase because of surface chemistry: deuterium is 
efficiently locked into grain-surface HDO. Efficient HD deple- 
tion, particularly in the very centers of the cores, means that the 
molecular ions previously thought to be tracers of the innermost 
regions of starless cores, H2D + and D2H + , will eventually dis- 
appear from the gas phase. It should be noted that in addition to 
these two ions, other deuterated species dependent on the abun- 
dance of H2D + , suchasDCO + andN2D + , show similar behavior. 
The timescale of HD depletion increases if chemical interaction 
between the gas and the dust precedes the starless core phase. In 
this case, part of the oxygen reservoir is locked in grain-surface 
H2O in the beginning of the starless core phase, which limits 
surface HDO production. 

While the present model predicts that HD ultimately depletes 
from the gas phase, our results do not imply that species such 
as H2D + and D2H + would be completely unusable as tracers of 
cold, dense gas. On the contrary, at later stages of core evolu- 
tion, the H2D + and D2H + abundances increase strongly towards 
the core edge (i.e. lower density) and the distributions become 
extended, which agree s with observation al evidenc e of extended 
distrib utions of H 2 D + dVastel et al.ll2006l in L 1544; |Pagani et all 
120091 in L183) and D7H + dParise et al.ll20lU in Oph H-MM1). 
In this sense, a high observed H2D + and/or D2H + abundance 
toward a starless or prestellar core could (loosely) constrain the 
core age, although the detailed behavior of the deuteration chem- 
istry is uncertain as it depends on various parameters, including 
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the initial H2(o/p) ratio. We note that the evaluation of the initial 
elemental abundances and the initial H2(o/p) ratio would require 
detailed modeling of the core condensation phase. Also, more 
laboratory work on the O + H2 surface reaction (and on tunnel- 
ing on grain surfaces in general) is required, as this reaction can 
significantly affect our results. 

A natural extension of this work is to perform a larger-scale 
study of deuterium chemistry including comparisons with obser- 
vations - this will be the subject of an upcoming paper. Finally, 
we note that the timescale of forming grain-surface HDO, criti- 
cal to the gas phase abundance of HD according to our results, 
depends on the surface abundances of OH and OD. As our model 
does not include a multilayer treatment of surface chemistry, 
it may overestimate the a bundance of radical s available on the 
surface (as pointed out bv lTaquet et al .1120121) . This could have 
implications on HD depletion, especially if HD depletion pro- 
ceeds through conversion into HDO* as in our models. While 
the present model does not allow for it, this issue should be in- 
vestigated. 
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fall into this category. For these reactions, we have assumed that 
only para H2 is formed (and that, accordingly, the rate coeffi- 
cient remains unchanged). Production of ortho-H2 requires an 
exothermicity of at least ~170 K (the energy difference between 
ortho and para states of H2). The exothermicity of each reaction 
can be calculated if the enthalpies of each species involved in the 
reaction are known - we have however not attempted such calcu- 
lations here, and proceeded on the assumption that the required 
exothermicity is not reached in general. 

A.3. Reactions involving H+ where H 2 is created 

The majority of the ion-molecule reactions of H^ are of the same 
type as reaction (O in the main text where H3 donates a proton 
and spin selection rules are applied. However, there are also re- 
actions such as 



H+ + MgH 



•Me 



H? + Ho 



(A.6) 



where multiple H2 molecules are formed. For these reactions, 
we assume the following separation rules: 



H+(o) + MgH -A Mg + + H 2 (o) + H 2 (p) , 



(A.7) 



Appendix A: Details of the ortho/para separation 

In this appendix we further discuss the ortho/para separation of 
the OSU reaction set. As mentioned in Sect.|2] there are several 
reaction types where branching ratios based on spin selection 
rules, such as in reaction (0, are not used; we list these special 
cases here. For each reaction type, we describe the assumptions 
behind the adopted branching ratios (if any) and give an example 
reaction. 



A.1. Charge transfer reactions 

We assume that in charge-transferring reactions such as 



kk A 



/<A 



H^+C 2 H-^C 2 H + +H2 



(A.1) 



the spin state of HJ (or H^ ) is conserved (we require that total 
nuclear spin in conserved in the reaction). The above reaction 
then separates into two reactions: 



H+(o) + C 2 H -H C 2 H + + H 2 (o) 
and 

H+(p) + C 2 H -H C 2 H + + H 2 (p) , 
with the same rate coefficient &ai ■ 



(A.2) 



(A.3) 



A.2. Reactions of species other than HJ or H+ where H 2 is 
created 

In these reactions, neither H} nor HJ appears as a reactant; thus, 
reactions such as 



H^(p) + MgH — * Mg + + H 2 (o) + H 2 (p) , 



and 



kka 



H 3 + (p) + MgH -^ Mg + + H 2 (p) + H 2 (p) . 



(A.8) 



(A.9) 



This logic is also applied to the (few) similar reactions that in- 
volve HJ and a molecule with several hydrogen atoms. In prin- 
ciple, reaction ( IA.6I 1 can also produce H2(o) + H2(o) which can 
be sh own b y calcula ting the branching rati os with the method of 
Oka (l2004t see also IWirstrom et aLll2012l) . However, branching 
ratios based solely on nuclear spin statistics might not be appro- 
priate for reactions with low exothermicies, and a simplifying 
assumption about the branching ratios is utilized. We note that 
while reactions such as ( IA.6t are not among the most impor- 
tant reactions controlling the spin states of H2 and H 3 in starless 
cores, this important issue warrants further investigation. 

A.4. Reactions involving H3 or Mi, with no H 2 produced 
These are reactions of the same type as the reaction 



H+ + 0- 



• OH + + H . 



(A.10) 



In these cases, ortho and para states are destroyed with equal 
rates, so that the above reaction separates into 



Hj(o) + O -A OH + + H , 
and 

Hj(p) + O -H OH + + H . 



(A.ll) 



(A. 12) 



NH + + H 2 -^ HNO + + H 2 

and 

C2H3 + H — > C2H2 + H2 



(A.4) a_$_ Other reactions 



(A.5) 



In addition to the reaction types presented above, there are also 
a number of reactions that do not fall into the general categories 
discussed above. Whenever this is the case, we apply a custom 
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logic unique to each reaction that is, as far as possible, consistent 
with the above general cases. For example, the reaction 



k ts 

HJ + H 2 S -A HS + + H + H 2 

separates into the reactions 

H+(o) + H 2 S — > HS + + H + H 2 (o/p) , 

and 

H+(p) + H 2 S — > HS + + H + H 2 (o/p) , 



(A.13) 



(A. 14) 



(A. 15) 



with 1 : 1 branching ratios for bo th reactions. For the former reac- 
tion, the method of lOkal d2004l) yields an H 2 ortho/para ratio of 
2: 1 (and 1 : 1 for the latter reaction), when H 2 S can exist in either 
ortho or para form. But since we do not know the spin state of 
H 2 S, we assume that it exists totally in para form and make the 
above simplifying assumption. 
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